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ABSTRACT 

Non-ideal MHD effects play an important role in the gas dynamics in protoplanetary disks (PPDs). 
This paper addresses its influence on the magnetorotational instability (MRI) and angular momen- 
tum transport in PPDs using the most up-to-date results from numerical simulations. We perform 
chemistry calculations using a complex reaction network with standard prescriptions for X-ray and 
cosmic-ray ionizations. We first show that no matter grains are included or not, the recombination 
time is at least one order of magnitude less than the orbital time within 5 disk scale heights, justifying 
the validity of local ionization equilibrium and strong coupling limit in PPDs. The full conductivity 
tensor at different disk radii and heights is evaluated, with the MRI active region determined by 
requiring that (1) the Ohmie Elsasser number A be greater than 1; (2) the ratio of gas to magnetic 
pressure /3 be greater than f3^i^{Am) as identified in the recent study by Bai & Stone (2011), where 
Am is the Elsasser number for ambipolar diffusion. With full flexibility as to the magnetic field 
strength, we provide a general framework for estimating the MRI-driven accretion rate M and the 
magnetic field strength in the MRI-active layer. We find that the MRI-active layer always exists at 
any disk radius as long as the magnetic field in PPDs is sufficiently weak. However, the optimistically 
predicted M in the inner disk (r = 1 — 10 AU) appears insufficient to account for the observed range of 
accretion rate in PPDs (around IO^^Mq yr~^) even in the grain-free calculation, and the presence of 
solar abundance sub-micron grains further reduces M by one to two orders of magnitude. Moreover, 
we find that the predicted M increases with radius in the inner disk where accretion is layered, which 
would lead to runaway mass accumulation if disk accretion is solely driven by the MRI. Our results 
suggest that stronger sources of ionization, and/or additional mechanisms such as magnetized wind 
are needed to explain the observed accretion rates in PPDs. In contrast, our predicted M is on the 
order of 1O~^M0 yr~^ in the outer disk, consistent with the observed accretion rates in transitional 
disks. 

Subject headings: accretion, accretion disks — instabilities — magnetohydrodynamies (MHD) — meth- 
ods: numerical — protoplanetary disks — turbulence 



1. INTRODUCTION 

Protoplanetary disks (PPDs) around pre-main- 
sequence stars form from the collapse of protostellar 
cores as a result o f angular momentum conservation 
(jAdams et ia,l.lll987D. With a typical lifetime of 1 - 10 
Myrs (e.g., lHillenbrand et al.lll998l: ISiciha-Aguilar et al.l 
|2006[ ). PPDs feed gas onto the central protostar, power 
an outflow and/or jet, and provide the raw materials for 
the formation of planetary systems. The structure, evo- 
lution and dispersal of PPDs are of crucial importance 
in understanding a wide range of physical problems 
es pecially in the a rea of planet formation (see review 
by lArmitagd 120111 ). As more and more exoplanets are 
discovered^, together with the advan cement of plane t 
formation theo ry (e.g., see the bo o k bv lArmitagd ()2010D . 
and also see IChiang &: YoudinI (|2010l ) for a review 
on planetesimal formation), understanding the gas 
dynamics in more detail in the PPDs becomes essential. 

One of the most important observational constraints 
relevant to the gas dynamics in PPDs is that PPDs 
are actively accreting. The accretion signature comes 
from the UV excess emission that veils the intrin- 



sic photospheric spectrum of a YSO, which is inter- 
preted as coming from the standing accretion shock 
formed at the s t ellar surface (jCalvet fc Gullbriii^ 119981 : 
iGullbring et al.l l2000l ). More commonly, the accre- 
tion rate M can be inferred from emission line pro- 
flies, in particular the H a line, based on the magne- 
tospheric accretion model (iMuzerolle et al.lll998l |2(301.1 . 
The inferred M f or Classical T-Tauri s t ars is about 
lQ-8±iM^ vr-i (JHartm ann et al.l 119981 : ICalvet et aH 
12004 IsTcilia-Aguilar et a l. 2005), and observations over 
a wide mass range of protostars reveal a correlation 
betw e en M and the pro t ostellar mass ([Muzerolle et al.l 
2005t iNatta et al.l I200I iHerczee fc Hillenbrand l2'008r 
Fang et al.l |2009() . Transitional disks, characterized by 
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optically thin inner holes whi ch may represent a later 
evolutionary stage of PPDs (jCa lvet et ajj I2002L 120051 : 
iHughes et ahl 120071: lEspaiflat et"al. 2007il , are also ob- 
served to be actively accreting. Although transitional 
objects are much rarer, the median of their accretion rate 
from currently available samples is a few times 10~^ M(^ 
yr-i (jNaiita et al.l 120071: Isicilia-Aguilar et al.ll2010D . A 
key question here is, what drives the rapid accretion in 
PPDs? 
The answer may depend on the evolutionary stages of 
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the PPDs, but it appears certain that turbulence plays 
a crucial role, and evidence of turbulence in PPDs has 
been re ported from sub-milh rneter interferometric obser- 
vations (jHughes et al.ll2011l ). In early phases, the accre- 
tion is likely to be primarily driven by the gravitational 
instability (GI). The GI leads to gravitoturbulence if the 
cooling rate is less than the orbital frequency, and trans- 
ports angular momentum out ward via non-axisymmetric 
spira l waves very efficiently (jGammid l2001t iRice et aLl 
|2005[) . For typical accretion rate in PPDs, and given the 
opacity dominated by dust grains, gravitoturbulence is 
likely to be present at i ntermediate disk radii between a 
few t e ns and ~ 100 AU (jVorobvov fc Basull2007t iRafikovl 
120091 : [Rice et al.ll2010t) . However, time-dependent cal- 
culations of the PPD evolution indicate that after the 
envelope infall stops (< 1 Myr), the disk is generall y 
not massive enough to sustain the GI ()Zhu et al.ll2010[ ). 
Other driving mechanisms at the inner disk as well as 
beyond the infall (embedded) phase are clearly needed. 

The most important mechanism for driving accretion 
in non-self-gravitating thin disks is believed to be the 
magn etorotational instability (MRI, iBalbus fc Hawlevl 
Il991| ). which generates turbulence and efficiently 
transports a.ngular momentum rad ially outward 
(iHawlev et alJ fl99l iStone et all [l996h . The MRI 
requires sufficient coupling between the gas and the 
magnetic field. However, PPDs are only weakly ionized. 
The main ionization sources such as cosmic rays and 
X-rays from the protostar only effectively ionize the sur- 
face layer of the disk, making the surface layers "active" 
to MRI driven turbulence, while the gas in the midplane 
remains poorly coupled to th e magnetic field and is 
termed 'dead" (|Gammielll996[ ). The layered accretion 
scenario has been studied in great detail via both numer- 
ical calculations in a fixed disk model with a chemical 
reaction network (jSano et al.l l2000t IFroman g et aQ 
[200l Hlmer fc NelsonI " IMifil: 
I2009D. and 
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ISemenov et al.l . — 

Bai fc Goodm anI 120091: [Turner fc Drake 



MHD s imulations (iFleming fc Stonell2003t [Turner et al.l 
20071: [iw TerfcSanol 120081 lllgner fc NelsonI 12001 



Oishi fc Ma c Low 200Q). These research works generally 
confirm that a "dead zone" is expected at the inner part 
of the PPDs (about 0.5 - 10 AU), though its radial and 
vertical extent depends on the ionization rate and the 
abundance of small (sub-micron) dust grains. 

Because of the low ionization level in PPDs, the 
gas dynamics is controlled by a number of non-ideal 
MHD effects, including Ohm ic resist i vity, Hall cfi'ect 
and ambipolar diffusion (see iBalbusI ([20111 ) for a re- 
view) . All calculations and simulations mentioned above 
considered only the effect of Ohmic resistivity. How- 
ever, Hall effect and ambipolar difi'usion (AD) also play 
an important role, and dominate Ohmic resistivity in 
the more ten uous and more strongly magneti z ed disk 
uppe r layers ([Salmeron fc Wardld I2005L [20081 : IWardld 
[2007() . The linear dispersion properties of the MRI in 
the Hall and AD regimes differ substantially from those 
in the Ohmic regime. In particular, the inclusion of 
the Hall effect makes the properties of th e MRI depend 
on the orientation of th e magnetic field (fWardld [19991 : 
IBalbus fc Teroueml[200l[ : iWardle fc Salmeronll201lin n 
the presence of both vertical and azimuthal field, MRI 



can grow at appreciable rate even i n the limit o f infinitely 
strong AD (|Kunz fc Balbudl2004l : lDeschll2004| ). 

The effects of Ohmic, Hall and AD on the non-linear 
evolution of the MRI have been studied separately by 
various numerical simulations. One of the key questions 
addressed by these simulations is that under which condi- 
tions can the MRI be sustained. It has been found that 
for Ohmic resistivity, MRI can be sustained when the 
vertic al Elsasser number A ^ = v\,/rioQ. is grea ter than 
unity ([Turner et al.l [20071: lllgner fc NelsonI |2008[ ) . where 
VAz is the vertical component of the Alfven velocity, rjo 
is the Ohmic resistivity, and J7 is the angular frequency 
of the disk. Numerical simulations including both the 
Ohmic resistivity and the Hall term were performed by 
ISano fc Stoni ([2002a[[bl ) . They found that the HaU effect 
does not strongly affect the saturation level of the MRI, 
and the A^ > 1 criterion for the MRI to be sustained still 
holds with the inclusion of the Hall term. Although their 
exploration of the Hah par ameter was not quite complete 
([Wardle fc Salmeronll201lh . these results are not surpris- 
ing given the fact that t he Hall term i s not dissipative. 

In our recent paper, [Bai fc Stond ([201 it ) (hereafter, 
BSll), we studied the effect of AD on the non-linear 
evolution of the MRI in the strong coupling limit, which 
applies when the ion inertia is negligible and the recom- 
bination time is much shorter than the orbital time. This 
assumption will be justified in this paper in the context of 
PPDs. The effect of AD is characterized by the Elsasser 
number based on AD. denoted by Am^ which describes 
the number of times a neutral molecule "collides" (i.e., 
exchanges much of its momentum) with the ions in an 
orbital time. The key result of BSll is summarized in 
their Figure 16: MRI can be sustained at any value of 
Am, provided that the magnetic field is sufficiently weak 
/3 > Pinin{Am), where /3 is the ratio of gas to magnetic 
pressure (see Equation (25) of BSll), and /3niin(^m) in- 
creases with decreasing Am. This relation provides an- 
other constraint on the sustainability of the MRI. 

In this paper, we aim at studying the location and ex- 
tent of the active regions in PPDs, and predicting the 
MRI-driven accretion rate in the most realistic man- 
ner by incorporating all the currently available numer- 
ical simulation results. We do so by solving a com- 
plex set of chemical reaction network established in 
[Bai fc GoodmanI ([2009t ) (hereafter BG09). A single pop- 
ulation of dust grains is also included in the network. 
Magnetic diffusion coefficients are calculated from the 
equilibrium abundance of charged species. A unique fea- 
ture in our treatment is that we have included the full 
dependence of magnetic diffusion coefficients (hence the 
Elsasser number) on the magnetic field strength, with the 
field strength constrained by the results from non-ideal 
MHD simulations. This allows us to predict the mag- 
netic field strength and the accretion rate in PPDs using 
the least amount of assu mptions. One closely related 
work is bv IWardld ([20071 ). who performed similar chem- 
istry calculations to obtain magnetic diffusivities of all 
non-ideal MHD effects with a simpler reaction network, 
but the extent of the active layer and accretion rate was 
not addressed in detail which may be partly due to the 
unavailability of numerical simu lation results. Another 
closely related work to ours is bv lPerez-Becker fc Chiang 
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(|2011af) . who were motivated by the accretion problem 
in transitional disks and the role of tiny grains. They 
have considered both Ohmic resistivity and AD, although 
their adopted criteria were more simplified and did not 
account for the role of magnetic field strength. 

This paper is structured as follows. We begin by re- 
viewing the derivation of various non-ideal MHD effects 
in Section [2] We describe our chemical reaction net- 
work and calculation of magnetic diffusivities in Section 
[31 where we also discuss our adopted criteria for the MRI- 
active layer. In Section |4] we present the results of our 
fiducial model calculation, where a framework for esti- 
mating the accretion rate and the magnetic field strength 
is provided. Using this framework, we study the depen- 
dence of accretion rate on various ionization and disk 
model parameters in Section [5] We summarize and con- 
clude in Section [5] 

2. OVERVIEW OF NON-IDEAL MHD EFFECTS 

Non-ideal MHD effects derive from the generalized 
Ohm's law. In the single-fluid framework for weakly ion- 
ized gas, fluid density p and velocity v specify the den- 
sity and velocity of the neutrals. Let charged species j 
has particle mass rrij, charge ZjC, number density Uj, 
and drift velocity relative to the neutrals Vj. Charge 
neutrality condition applies for non-relativistic MHD: 
^ ■ UjZj = (note that Zj can be either positive or nega- 
tive). Let E' be the electric field in the frame co-moving 
with the neutrals, while for non-relativistic MHD, the 
magnetic field B is the same in all frames. In this co- 
moving frame, the equation of motion for charged species 
(whose inertia is negligible) is set by the balance between 
the Lorentz force and the neutral drag, given by 



Z,e{E' - 



X -B) = IjPmjVj 



(1) 



where 7j = {av) j / {m + mj) with {av)j being the rate co- 
efficient for momentum transfer between charged species 
j and the neutrals, and m is the averaged particle mass 
of the neutrals. 

The relative importance between the Lorentz force and 
the neutral drag is characterized by the ratio between the 
gyrofrequency and the momentum exchange rate 
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Charged species j is strongly coupled with neutrals if 
\/3j\ <^ 1, and is strongly tied to magnetic fields when 
1/3,1 »1. 

Since the current density is given by J = e^ ■ rijZjVj. 
The generalized Ohm's law can be obtained by inverting 
equation ([T]) to express Vj as a function of E' . The result 
is 



J ^aoE\,+aHB X E'^ + apE'^ 



(3) 



where subscripts || and j. denote vector components par- 
allel and perpendicular to the magnetic field B, and 
denotes unit vector. The Ohmic, Hall and Pcderscn con- 



ductivities are (jWardlell2007l ) 



ao = -^^njZjPj , 
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respectively. Note that the Hall conductivity depends on 
the sign of Zj , while ctq and ap are always positive since 
Pj has the same sign as Zj. 

The Ohm's law (|3]) can be inverted to give the elec- 
tric field using current densities, which then leads to the 
induction equation modified by non-ideal MHD terms 

dB 47r 

— - = V X (i; X B) V X [?7oJ + ??h(J x B) + 77aJ_l] , 

ot c 

(5) 

where 
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are the desired Ohmic, Hall and ambipolar diffusivities 
as in equation ([5]), determined by the microphysics of 
ion-neutral and electron-neutral collisions, and a±_ = 
\/crjj + ap . Note that only rjo is independent of B. 
The absolute value of these diffusion coefficients deter- 
mines the relative importance of the Ohmic, Hall and 
AD terms. 

The most commonly used magnetic diffusivities are ob- 
tained by assuming that electrons and positively charged 
ions are the only charge carriers, with all ions having 
the same Hall parameter. In this case, one can express 
the conductivities in terms of /3e and (3i for electrons 
and ions respectively, and since |/3e| ^ /3i, one finds 
(|Salmeron fcWardld 120031) 



dB 



Vx{vxB)-Vx 
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(7) 

Physically, because the electrons are the most mobile 
species in the gas, the magnetic field is effectively car- 
ried by the electrons. Correspondingly, the Hall and AD 
terms originate from electron-ion drift and ion-neutral 
drift respectively. In dense regions with weak magnetic 
field, both electrons and ions are well coupled to the 
neutrals (1 ^ |/3e| ^ Pi)-, and Ohmic resistivity dom- 
inates. In tenuous regions with strong magnetic field, 
both electrons and ions are tied to the magnetic field 
(|/3e| ^ Pi ^ 1), and AD dominates due to the ion- 
neutral drift. The Hall dominated regime due to the 
electron-ion drift lies in between, where electrons are tied 
to the magnetic field while the ions are coupled to the 
neutrals {\pe\ ^ 1 ^ Pi)- The above formula no longer 
holds when a substantial fraction of charged particles are 
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grains (|Baill201ll) . although it is widely used for its sim- 
plicity and physical clarity. In this paper, we shall adopt 
the general expression for magnetic difFusivities ^. 

The energy dissipation rate associated with the non- 
ideal MHD effects is given by 



E = -En 

c 






(8) 



where £„ = (47r/c)[?7oJ + i]h{J x B) -f tjaJ i\ is the 
electromotive force associated with the non-ideal MHD 
terms. We see that the Ohmic resistivity dissipates the 
total current (leading to magnetic reconnection) , while 
AD damps the perpendicular component of the current 
(via ion- neutral drag). On the other hand, the Hall effect 
is not dissipative. It describes the magnetic diffusion due 
to drift motion between charged carriers without break- 
ing magnetic field lines, and is also present in fully ionized 
plasma. 

3. CALCULATIONS OF NON-IDEAL MHD EFFECTS IN 
PPDS 

As discussed in Section [5J full assessment of the non- 
ideal MHD effects requires knowledge of the number den- 
sities of all charged species in PPDs. In this section, we 
describe our chemistry calculations to infer magnetic dif- 
fusivities in PPDs and provide a quantitative criterion for 
judging whether the MRI can operate or not under the 
given diffusivities based on results from numerical sim- 
ulations. Most of the chemistry calculation procedures 
are adopted from BG09. 

3.1. Disk Model 

Fiducially, we take the minimum-mass solar neb ula 
(MMSN) model (IWeidenschillindUOTTt iHavashilfTosl as 
our disk model, which is simply constructed by smearing 
the required mass for forming the solar system planets 
into a smooth distribution. It represents the minimum 
amount of disk mass required for forming the solar sys- 
tem planets. The surface mass density and disk temper- 
ature are given by 
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(9) 



where tau is the disk radius measured in AU, and the 
disk is treated as vertically isothermal. By default, we 
assume the mass of the protostar to be A'/* = IMq, with 
the Keplerian frequency 51 = \J GM^jr'^ . The mean 
molecular weight of the neutrals is taken to be /i„ = 
2.34771//, from which the sound speed Cg = y^ kT / finirip 
and disk scale height H = Cs/il can be found. 

In additio n, we also co nsider the solar nebula model 
proposed by iDeschI ()2007() , which takes into account the 
recent advances in the pla net formation theory (in par- 
ticular, the "Nice" model, iTsiganis et al.l [2005[) . and is 
much more massive than the MMSN. The surface den- 
sity and temperature profiles are given by 

I]„ = 5 X 10V7,V^ g cm-2 

-0.43 .. (10) 



T = 150r 



AU 



K, 



where the temperature profile is estimated from 
iChiang fc GoldreichI (Il997l) . 



Submillimeter interferometric observations of the 
^ 1 Myr ol d Oph i uchus star forming regions by 
[Andrews etldl (|2009l |2010[ ) have revealed the density 
and temperature profiles in the outer regions (> 10 AU, 
due to limited spatial resolution) for a sample of PPDs. 
The surface density for the majority of the disks appears 
to match the MMSN value well at 10 - 20 AU. Although 
the fitted density profile in the outer disk is shallower 
than the MMSN and the Desch's model (with the me- 
dian slope of —0.9 rather than —1.5 or —2.2), the surface 
density profile of the inner disk is not well constrained 
by observations, and both MMSN and the Desch's model 
may be viable choices. Furthermore, global calculations 
of PPD evolution do indicate much higher surface mass 
densities in the inner disk than direct c ontinuation of th e 
observed density profiles to small radii ()Zhu et al.ll2010|) . 

3.2. Ionization Sources 

The PPDs are generally too cold for thermal ionization 
to take place except in the innermost regions (< 1 AU, 
iFromang et al.|[2002[ ). We are interested in regions with 
r > 1 AU and consider the following three non-thermal 
ionization sources. 

First, the X-ray ionization from the protostar. Most T- 
Tauri stars produce stro ng X-ray emission due to corona 
activities (see review by iFeigelson et al.l |2007() . The X- 
ray fluxes are generally variable, with larg e X-ray flares 
recurring on time scale of a few weeks (iStelzer et al.l 
|2007() . The X-ray emission during the flares is harder 
than that in the quiescent state. The time averaged X- 
ray luminosity is roughly proportional to stellar mass, 
and is about 10^^ to 10^^ erg s~^ for solar mass stars 
(jPreibisch et al.ll2005l: [Gudel et al.ll200 7^, with typical X- 
ray temperature ranging from 1-8 keV (Wolk et al. 2005). 
We adopt the X-ray temperature Tx = 5 keV, and X-ray 
luminosity Lx = 10'^° erg s~^ as our standard model pa- 
rameters. We tak e the io nization rate (^3?) calculated by 
llgea fc Glassgoldl (|1999[ ). which takes into account both 
absorption and scattering of X-ray photons. In practice, 
we use the fitting formula given by BG09 (see their Equa- 
tion (21)) 



toff 
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(11) 



where Lx,30 = Lx/10^°erg 8"^ TVi = 3.0 x lO^i cm'^, 
7V2 = 1.0 X 10^* cm~^, Nh denotes the column number 
density of hydrogen nucleus from the point of interest 
to one side of the disk surface, while the "bot." symbol 
represents the dual terms with Nh being the hydrogen 
column number density to the other side of the disk sur- 
face. The column number densities TVi and iV2 roughly 
correspond to column mass density of 1.2 x 10^^ g cm~^ 
and 3.9 g cm~^ respectively. 
Secondl y, the cosmic-ray (CR) ioniza tion with ioniza- 



tion rate (jUmebavashi fc Nakanci|[l98lh ^ 
C^l = 1.0 X 10"^^ exp (-S/96g cm^^) g-^ + bot. (12) 

^ See lUmebavashi &: Nakanol I I2009I ) for a more refined formula. 
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The CR flux is highly uncertain because on the one hand, 
the flux may be much higher if a supernova explosion oc- 
curs in the vicinity of the protostar, and observations of 
the CR flux toward the diffuse cloud C Persei indicate en - 
hanced ionization rate of lO^^*^ s"^ (|McCan et alJl2003[ ): 
but on the one hand, the CR flux may be substantially 
shielded by the stellar wind. 

Third, the radioactive decay, primarily the decay of 
short lived ^^Al, produces ioniz ation rate of 3.7 x 10~^ ^ 
s-i with half-life 0.717 Myr (jTurner fc Drake! I2009D . 
Here we adopt the ionization rate of 10""'^^ s~^ as ap- 
propriate for disk ages of around 3 Myr. The radioactive 
decay is generally too weak to provide sufficient ioniza- 
tion, but it prevents the midplane of the inner disk from 
being completely neutral. 

Other possible ionization source s such as energetic pro- 
tons from disk and stellar corona (jTurner fc Drakell2009t ) 
are ignored for simplicity. Although the resulting ioniza- 
tion rate may exceed X-r ay and cosmic-ray ionization b y 
a factor of as large as 40 (jPerez-Becker fc Chiand[2011a[ ) , 
their fluxes are highly variable and uncertain. In our cal- 
culations we also consider Lx = 10'^^ erg s~^ whose ion- 
ization rate is likely to overwhelm this effect by at least 
an order of mag nitude. 

Very recently, iPerez-Becker fc Chiand (|2011bf ) pointed 
out another potentially important ionization source: the 
far ultraviolet (FUV) radiation from the protostar. FUV 
photons are unattenuated by the hydrogen column, and 
efficiently ionize tracer species of heavy elements such 
as C and S with penetration depth of up to O.lg cm~^. 
FUV ionization is not included in our calculation, while 
its relative importance will be discussed in Sections 14.31 
and 



mass distribution /(a) (with J f{a)da = Z = 0.01), one 
may conveniently consider 



3.3. Grain Size and Abundance 

The abundance and size distribution of grains in PPDs 
are crucial for disk chemistry. Observationally, they are 
constra ined by mode ling the spect ra energy distribution 
(SEP) (IChiang fc Go ldreich 1993 ID'Alessio et "all[l998l: 
iDuUemond fc DominikI 120041 ). Although the parameters 
are very degenerate, there ha ve been evidence of g rain 
growth to above micron size (jD'Alessio et al.ll2001|) . as 
well as dust settlin g (to the midplane , which exhibits 
as grain depletion IChiang et al.l 120011: JD'Alessio et al.l 
120061: iFurlan et all 120061: IWatson et al.l I2009D . Results 
from mid-infrared spectro scopy from both T-Taur i stars 
and Herbig Ae/Be sta rs (|van Boekel et al.l l2003l 120051: 
iPrzvgodda et al.l |2003[ ) also revealed the presence of 
micron-sized grains in a substantial fraction of PPDs. 
In our chemistry calculation, we adopt a simplified pre- 
scription of single-sized, well-mixed grains, and consider 
grain sizes of . 1 /im and 1 ^m . The internal density of the 
grains is taken to be 3 g cm~'^. Although a size distribu- 
tion of grains and some level of grain settling would be 
more realistic, our treatment is sufficient to demonstrate 
the basic physics. Moreover, BG09 considered two popu- 
lations of grains and found that the two populations be- 
have essentially independently. They further found that 
for a range of grain sizes considered (a = 0.01/im to 
l^m), the controlling parameter lies somewhere between 
the total grain surface area, and the grain abundance 
weighted by linear size. For a population grains with 
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0.01 



(13) 



as a measure of the chemical significance of dust grains. 
Our choice of a = 0.1/im and a = Ifim correspond to 
S* = 32 and 5 = 1 respectively. For a continuous size 
distribution of grains, the situation can be more compli- 
cated, but one may calculate the resulting S and use our 
results as a guide. 

Another potentially important ingredient of dust 
grains is the polycyclic aromatic hydrocarbon (PAH), 
which represents the smallest end of grain size dis- 
tribution. More than 50% of the Herbig Ae/Be 

disks have been detec ted to have PAH emission 

(jAcke fc van den Anckerl |2004[) . while the frac- 
tion i n T-Tauri sta r s is much smaller (jGeers et al.l 
120061 lOliveira et all |2010D . The importance of 

PAHs on the disk co n ductivi ty has been raised by 
iPerez-Becker fc Chiang! (j2011a| ) recently. They argue 
that PAHs in T-Tauri disks may be equally abundant as 
in Herbig Ae/Be disks, but they are undetected because 
T-Tauri stars are much less luminous in ultraviolet 
(UV) radiation to excite fluorescence emission from the 
PAHs. The abundance of PAHs in PPDs is not well 
determined, but the small size of PAHs means that they 
can have very large abundance (which rapidly recombine 
electrons) without contributing much to the total dust 
mass, posing a potential threat to the MRI. In this 
paper, we do not include PAHs in our chemical network 
by assuming they are ins ufficiently abundant. However, 
in our companion paper ()Bail!201ll ). we point out that 
due to their own conductivity, charged PAHs can even 
facilitate the MRI by suppressing the Hall effect and 
ambipolar diffusion under certain situations. 

3.4. Chemical Reaction Network 

We adopt the comp lex chemical reaction n etwork de- 
scribed in full detafl in !llgner fc Nelson! (!2006[ ) and BG09 
(see their Section 3). It contains nine elements, 174 gas- 
phase species and 2083 gas-phase chemical reactions ^. 
The rate coefficients are taken from the UMIST database 
(jWoodall et al.l[2007! ) (except for the ionization reactions, 
whose rate is given in Table 1 of BG09). We further in- 
clude a single population of grains as described in the 
previous subsection. Grain related reactions include col- 
lisional charging with electrons and ions, adsorption and 
desorption of neutrals, and grain collisions. The max- 
imum grain charge is taken to be ±10 and ±30 ele- 
mentary charges for O.lfim and 1/im grains respectively. 
This leads to substantial increase of total number of 
species (to 266 for a = O.lfim) and total number of re- 
actions (to 4513 for a = 0.1/im and more than 9000 for 
a = 1/xm) compared with BG. We have made a correction 
to !llgner fc Nelson! ()2006l) and BG09 in the calculation of 
the electron sticking coefficient, which is described in Ap- 
pendix |X] and discussed in Section UTTJ The calculation 

^ We correct the number of gas-phase reactions in BG09 where 
it was stated to be 2113. 
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of all other reaction rate coefficients remain unchanged 
from BG09. 

We initialize the chemical reaction network from single- 
elem ent species, whos e abun dance is provided in Table 
6 of lllgner fc NelsonI ()2006l ). For the two metal ele- 
ments Mg and Fe, we take their abundance relative to 
hydrogen nucleus to be 1.0 x 10~* and 2.5 x 10~^ re- 
spectively (same as most calculations in BG09). The 
grains are assumed to be uniformly mixed in the disk 
with 1% in mass. The large set of stiff ordinary differ- 
ential equations are evolved by the fourth-order implicit 
Kaps-Rent rop integrator with adaptive time stepping de- 
scribed in iPress et al.l (|1992[ ). Conservation of charge 
and elemental abundance is enforced at each step of evo- 
lution. The chemical network is evolved for 10^ years 
when quasi-cquilibrium has reached. 

3.5. Calculation of Magnetic Diffusivities 

To formally obtain the conductivities one needs to 
know the momentum transfer rate coefficients < av >. 
For ion-neutral collisions, the neutral atom is induced 
with an electrostatic dipole moment as it approaches the 
ions, the resulting rate coefficient is approximately in- 
dependent of tempera ture, and is in versely proportional 
to the reduced mass (jDraind (|201lD . see Table 2.1 and 
Equation (2.34)): 



< a« >,= 2.0 X 10-3 f^^ 



1/2 



T —1 

cm s 



(14) 



where fi = rtii^inlimi 4- /i„) is the reduced mass in a 
typical ion-neutral collision. 

For electron-neutral colli sions, we adopt the approx- 
imate fitting formula from iDraine et al.l ()1983[ ) , which 
applies at T > lOOK. The dependence of the rate coeffi- 
cient on T becomes shallower at lower temperatures due 
to the polarization effects, and we adopt 



< av >p= 8.3 X 10 X max 



1 



T 



lOOA' 



1/2-1 



"\ —1 

cm s 



(15) 
as an approximation. 

For collisions between neutrals and charged grains, the 
rate coefficient follows equation ()14|) for sufficiently small 
grains, while the collision cross section becomes geomet- 
ric for large grains. Therefore, we have 



< av >qr= max 



1.3 X 10" 



1.6 X 10" 



_i f a \ / T \ 



1/2. 
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With these collision rate coefficients, it is then straight- 
forward to evaluate the magnetic diffusion coefficients 
ivo, Vh and tja) from equations ([2]), (|4]) and dH). 

3.6. Recombination Time 

The recombination time ircb is an important quantity 
for studying the gas dynamics in PPDs. If trcb is much 
shorter than the dynamical time (r2~^) , as is requ ired in 
the so called "strong coupling" limit (jShul 119911 ). local 
ionization equilibrium would be a good approximation 



and the magnetic diffusion coefficients can be directly 
evaluated from the ionization rate and local thermody- 
namic quantities such as density and temperature. This 
will simplify numerical calculations considerably because 
a single- fluid approach is sufficient. In the opposite limit, 
if ^rcb is much longer than the dynamical time, a more 
appropriate approach would be the multi-fluid method. 
The recombination time ircb is not a well-defined quan- 
tity when there are multiple species of ions and when 
grains are present. Here we propose an effective recom- 
bination time t°^ by noticing the fact that resistivity 
scales linearly with the abundance of ionized species (es- 
pecially free electrons, see Section [2]). After evolving the 
reaction network to chemical equilibrium, we turn off the 
ionization sources and let the system relax (recombine) 
for a short period of time (^1 year). We measure the 
rate of change of the Ohmic resistivity, and the effective 
recombination time is defined as 

,off ^ Vo /-|7^ 

This definition captures the contribution from all recom- 
bination channels and is sensitive to the most rapid re- 
combination processes. It naturally generalizes the re- 
combination time deflned for single ion-electron recom- 
bination tj^lj = ne/{dne/dt) ex l/rie-^ 

3.7. Criteria for Sustaining the MRI 

In weakly ionized disks, non-ideal MHD terms dom- 
inate the inductive term if the Elsasser number is less 
than 1, where the Elsasser number based on the Ohmic, 
Hall and AD terms are defined as 



A: 



Am 



Vi 



rjofl 

-'a 



vz 



rjH^ 

A. 






(18) 



Here va = ^B^ /Anp is the Alfven velocity, and the ap- 
proximate equality in the definition of x and Am holds 
in the grain-free case. The Hall frequency LOh = eBn^/ pc 
is the cutoff frequency of the left polarized Alfven waves, 
which can be rewritten to ivu = {nf./n){mi/fin)'-^ci with 
uJci being the ion cyclotron frequency, n being the num- 
ber density of the neutrals, rui being the ion mass. 

Non-ideal MHD effects change the linear properties of 
the MRI substantially when any of the above Elsasser 



numbers falls below 1 (iBlaes fc BalbusI 



1991 IJml [19961: 



Kunz fc BalbusI 



Wardlel[T99l IBalbus fc TerauemI [2001 
20041 ). Below we summarize numerical studies on the 
non-linear evolution of the MRI in these non-ideal MHD 
regimes and provide our criteria on the strength and sus- 
tainability of the MRI turbulence, which are crucial to 
this work. 

For the Ohmic resistivity, vertically stratified shear- 
ing box simulations have identified that the border be- 
tween the MRI active and MRI inactive regions is well 

^ One can also define the efli'ectivc recombination time based on 
Hall and ambipolar diffusivities, which gives similar numbers but 
may have weak dependence on the magnetic field strength. 
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described by A = 1 (jllgner fc Nelson! |2008[ ). or A^ = 1 
([Turner et al.ll2007h . where the vertical Elsasser number 
Az = v\^/rio^, with vaz being the vertical component of 
the Alfven velocity. The former criterion gives slightly 
thicker active layer since A is generally larger than A^ 
by a factor of 3 — 30. The difference between these two 
criteria can be accommodated by noticing the the fact 
that transition from the MRI active to MRI inactive re- 
gions (i.e., the "undead zone") is relatively smooth, with 
the extent of the tra nsition region to be about Q.5H 
(jTurner fc Sancill2008[ ). Therefore, in our calculations, 
we will simply adopt A > 1 as the criterion for the ac- 
tive layer in the Ohmic dominated regime. It provides an 
optimistic estimate of the lower boundary of the active 
layer. 

For the Hall effect on the M RI, the only ex i sting n on- 
hnear study is performed by ISano fc Stond ()2002al lbl). 
They are motivated by whether the suppression of the 
MRI by Ohmic resistivity is affected by the Hall effect 
and performed simulations including both the Ohmic and 
Hall terms. They found that the saturation level of the 
MRI is not affected by the Hall effect by much. In par- 
ticular, the condition for sustained MRI turbulence is 
still controlled by the Ohmic Elsasser number and ap- 
pears to be independent of the Hall effect. However, the 
range of the Hall Elsasser number x (or X = 2/x in their 
notation) studied by Sano fc Stone is relatively narrow. 
Whether their conclusion can be generalized to the Hall 
dominated regime is still an open question. In this pa- 
per, we tentatively adopt Sano fc Stone's conclusion and 
ignore the Hall effect on the sustainability of the MRI, 
but we also discuss the applicability and potential caveat 
about this simplification in Section [4.21 

The effect of AD on t he non-linear evo l ution of the 
MRI has been studied by iHawlev fc Stond (|1998l ) using 
a two-fluid approach, applicable when the ionization frac- 
tion is large and the recombination time is long, and by 
BSll when the ion inertia is negligible and when ij-cb is 
much shorter than the orbital period (i.e., the strong cou- 
pling limit). The behaviors of the MRI in the two limits 
differ significantly. As wc will show in Section 14.21 the 
strong coupling limit applies almost everywhere in typ- 
ical PPDs. In this regime, BSll showed that the MRI 
can be self-sustained for any value of Am as long as the 
magnetic field is sufficiently weak. At a given Am, the 
maximum magnetic field strength is given by 



^n 
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V AtoO-S 



1 



1/2 



(19) 



where the plasma (3 = Pgas/Pmag is the ratio of gas to 
magnetic pressure, with Pmag = B'^ /8tt. Here (3 is not to 
be confused by the Hall parameter for charged particles 
/3j, which is used only in Section [2] 

In sum, our adopted criteria for sustained MRI turbu- 
lence is 

A > 1 , and /3 > /3,„i„(ATO) . (20) 

Besides the potential uncertainty from ignoring the Hall 
effect, another uncertainty in our adopted criteria is 
that they are mostly based on unstratificd simulations 
with vertical box size fixed at H (e.g., simulations by 
Sano & Stone and BSll). Although in the Ohmic 



regime unstratificd an d stratified simulations tend to 
yield similar crite rion JSano et al.l 119981 : iFleming et al.l 
I2000t [Turner et al.ll2007f ). it is vet to be explored whether 
the same situation holds for the case of ambipolar diffu- 
sion. In addition, above the MRI active layer, magnetic 
dissipation may generate a hot disk corona, which on the 
one hand, possesses some str ess (although much sm aller 
than that in the active layer. [Miller fc Stondl2000t ). and 
on the other hand, increases /3 in the upper disk. De- 
spite all these complications, our approach serves as the 
first step towards more realistic criteria, and moreover, 
it is sufficiently simple for illustrating our new method 
for estimating the location of the active layer and the 
accretion rate in Section [ 



3.8. Accretion Rate and Required Field Strength 

In the active layer, the MRI generates turbulent stress 
T,.0 which transports angular momentum outward, and 
its s trength is usually charact erized by the a parame- 
ter ([Shakura fc Sunvaev[ 119731 ). defined as T^^ = apc^. 
When MRI is self-sustained, there is a tight correlation 
between a and the time and volume averaged magnetic 
energy, which is ag ain characterized by the plasma (3 
([Hawlev et al.l[l995[ . BSll) 



a ; 



1 

2^ 



(21) 



Note that this relation holds only in the MRI-active layer. 
If accretion in PPDs is solely driven by the MRI in the 
active layer, a simple relation can be derived connect- 
ing the accretion rate M to the magnetic field strength 
(BG09). For steady state accretion, conservation of an- 
gular momentum demands 



Mnr^ = 2TTr^ 



dzT,, 



r0 



27rr^ 



dzTrc 



(22) 



active 



where the last integral is performed across the ac- 
tive layer. Although the dead zone has also been 
shown to transport angular momentum by non- 
axisymmetri c density waves launc h ed from the ac- 
tive layer (iFleming fc Stond l2003t lOishi fc Mac Low! 
I2009r). which is a generic proces s in shear turbulence 
( Heinemann fc Papaloizoull2009al [bl). its contribution is 
only a small fraction of that from the active layer, and 
can be safely ignored. From equation (|2ip . we have 
Trif, = aPgas ~ Pmag/2 in the active layer, and the above 
equation turns to 



M 



active °^' 



(23) 



This is a very useful formula for accretion rate estimation 
and is quite general for MRI driven angular momentum 
transport in both ideal MHD or non-ideal MHD regimes. 
It will be used extensively in Sections 14.31 and Section [SJ 
In turn, for MRI driven accretion, one can estimate the 
strength of the magnetic field given the accretion rate. 
Let the thickness of the active layer be denoted by ha for 
each side of the disk, the integral over the vertical height 
can be replaced by a factor of 2/iq, which leads to 



(B' 



AMn/ha 



(24) 
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where the bracket (•) means vertical averaging. To obtain 
a more quantitative estimate of the field strength applica- 
ble to PPDs, we consider the MMSN aromid a lAf© pro- 
tostar. The thickness of the active layer should generally 
be on the order of the disk scale height ha ^ H = Cs/^, 
and we obtain 



B) w 2jMn/H ; 



lomI/^a"/^ G 



(25) 



where M_8 = M/IQ-^Mq yr'^ We note that in ob- 
taining the above relation, only the temperature profile 
(to estimate the disk scale height) of the MMSN model 
is used (which is more reliable than the surface density 
profile) . This relation states that for MRI driven angular 
momentum transport, strong magnetic field is needed for 
fast accretion. 

For typical accretion rate of 10~^ Mq yr~^, the implied 
magnetic field strength is strong, and is in fact close to 
the equipartition strength for the active layer. Assuming 
the column density of the active layer Eq to be 10 g 
cm~^ (as comparable to the penetration depth of the 
X-ray ionization), the equipartition field strength is 
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where Si = Sa/10 g cm~^. The simple analysis here 
shows that to the limiting factor for MRI to drive ac- 
cretion in PPDs is not only the level of ionization, but 
also the strength of the magnetic field. A quantitative 
manifestation of this effect will given in Section 14.31 

4. ACTIVE LAYER IN THE FIDUCIAL MODEL 

We refer our fiducial model to the MMSN disk with 
X-ray luminosity of the protostar being lO'^^erg s~^ and 
with cosmic-ray ionizations. Variations to the fiducial 
model are discussed in the next section. Other parame- 
ters are fixed at values specified in the previous section 
in all calculations. 

Using the fiducial model, we run four calculations at 
two disk radii lAU and lOAU, with and without grains 
(0.1/im). In each calculation, we evolve the chemical net- 
work at fixed radius and scan from the disk midplane up 
to 5 disk scale heights. At each point, we extract the 
number density of all species from the end of the evolu- 
tion (10^ years) and calculate the magnetic diffusivities 
as a function of the magnetic field strength. Various 
aspects of the results are discussed in the following sub- 
sections. 

4.1. Chemistry 

In Figure [TJ we show the vertical density profile of 
various chemical species normalized to the number den- 
sity of the hydrogen nuclei for calculations at 1 AU. 
The most important quantity is the ionization fraction 
Xe = Ue/riH, as plotted in red, which largely determines 
the strength of the non- ideal MHD effects (Section [2]). 
In addition, we plot the profile of the ionization rate. 
This figure is to be compared with Figures 3 and 6 in 
IWardld (|2007( ). It is clear that the ionization fraction is 
extremely small (^ 10~^ in general), hence the inertia 
of the charged particles is negligible compared with the 
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Fig. 1. — Fractional abundance of electrons (red), metal ions 
(blue) and other ions (cyan) relative to hydrogen nuclei as a func- 
tion of height {z) above the midplane in our fiducial model at lAU. 
Also shown is the ionization rate (bold dashed) as a function of z, 
divided into three segments whore the dominant sources of ioniza- 
tion are labeled. Upper panel: calculation without grains. Lower 
panel: calculation with well mixed O.l/^m grains with 1% in mass. 
Abundance of positively charged (green), neutral (dark green) and 
negatively charged (magenta) grains are labeled by their elemental 
charge Z. 



inertia of the neutrals, justifying the first requirement of 
the strong coupling limit. 

The main driving force of chemical evolution is the ion- 
ization reactions. In the Figure, there are several "steps" 
in the vertical profile of the ionization rate that are as- 
sociated with transitions to different ionization regimes, 
as described in Section [X21 These "steps" also make the 
vertical profile of the ionization fraction x^ exhibit simi- 
lar features. At the uppermost layer, the ionization rate 
is the largest and is dominated by direct X-ray ionization 
from the protostar, with a very small column density of 
about 0.01 g cm^^. Slightly deeper down, the ioniza- 
tion is still dominated by the X-rays, but is mainly due 
to the Compton scattered X-ray photons from the upper 
layer, with penetration depth of about 4 g cm~^. Further 
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deeper, cosmic-ray ionization with penetration depth of 
about 100 g cm~^ takes over, while around the midplane, 
radioactive decay dominates. Comparing the result to 
Figure [2] shown in the next subsection reveals that the 
ionization rate provided by radioactive decay is generally 
too small to produce sufRcient ionization in PPDs, and 
can be safely ignored for the purpose of estimating the 
extent of the active layer and dead zone. 

In the grain-free calculation, we see that the ionization 
fraction primarily overlaps with the metal abundance at 
z < 4i/, indicating that the metals are the dominant 
electron do nor, which has been shown in many previous 
works (e.g., iFromang et "all 120021 : lllgner fc NelsonI I2006L 
BGG9). Above Mi, essentially all the metal atoms are 
ionized, and the main electron donor is taken over by 
other ions. The ionization fractio n from our calculation 
differs from that in lWardlg (|2007t ) mainly because we use 
a more complex (and presumably more realistic) chemi- 
cal network. 

The inclusion of well-mixed dust grains with a — 0.1/im 
dramatically reduces the ionization fraction. At the mid- 
plane, Xe is reduced by 5 orders of magnitude as com- 
pared with the grain-free case. The reduction factor is 
still significant but smaller in the upper layers up to 
z > hH . With grains, the role for metals as the main 
electron donor is suppressed because the recombination 
of metal ions is facilitated by grains, consistent with 
iWardld ll200l . 

The reduction of electron density by grains has an- 
other consequence: when the electron abundance falls 
substantially below the grain abundance, the ions and 
grains take over from the electrons to play a decisive 
role on the conductivity. In the chemistry calculation 
shown in the bottom panel of Figure [H this corresponds 
to regions close to the midplane with \z\ < 2.5H. The 
grain-free formula ([7|) for magnetic diffusivities no longer 
holds in this regime: the resistivity r]o is smaller than rje 
due to contribution from ions and grains. In addition. 
Am is no longer independent of magnetic field strength 
(as can be traced from Figure [3|) , and becomes larger in 
stronger field. This effect and its significance is discussed 
in full detail inHil (|2011[) . 

One differ ence between our calculation and the cal- 
culations by IWardld (|2007f ) is that we have considered 
the electron sticking probability Se- This probability 
Se was simply taken to be 1 in their calculation. In 
iPerez-Becker fc Chiang! ([20lI3), Se wa s fixed at 0.1 fo r 
PAHs, and 1 for normal grains, while in lOkuzumil (j2009( ). 
Se was fixed at 0.3 for all grain sizes. However, because 
electron is much lighter than the grain surface atoms, en- 
ergy transfer by inelastic collisions with grains is ineffi- 
cient and the electron may have a large chance to escape. 
The derivation of electron sticking c oefficient wi t h neu - 
trals grains has been performed by iNishi et all (|1991h , 
who showed that Se decreases with increasing temper- 
ature, from about 1 at zero temperature to about 1-2 
orders of magnitude le ss than unity at a few hundred K. 
lllgner fc NelsonI (|2006f ) and BG09 adopted this formula 
in their chemistry calculations, but they erroneously used 
the same formula for both neutral and charged grains. A 
generalized derivation of the electron sticking coefficient 
to include grain charges is given in Appendix [Aj where 



we find that the sticking coefficient for more positively 
charged grains is progressively smaller than that for nega- 
tively charged grains. This is mainly because the electron 
is accelerated more as it approaches the more positively 
charged grains, making it more difficult to get rid of the 
excess energy for adsorption. 

The inclusion of the sticking coefficient in grain- 
electron collisions reduces the electron recombination 
rate with grains, leading to less dramatic reduction of 
electron abundance and conductivity. Our test runs in- 
dicate that that the ionization fraction calculated with 
and without including the electron sticking probability 
can differ by up to a factor of 10 in certain regions. The 
electron sticking probability also affects the grain charge 
distribution. For example, in Figure [U the mean grain 
charge in the disk upper layers from our calculation is 
ab out —2. 5 elem entary charge rather than around —8 
in IWardQ (|2007f )'s calculation. Our new result on the 
charge dependence of Se implies that grains tend to be 
(slightly) more positively charged, although this is a rel- 
atively weak effect and is more relevant for sub-micron 
grains. 

4.2. Magnetic Diffusivities and Recombination Time 

The abundance of all charged species from the chem- 
istry calculation is used to evaluate the magnetic dif- 
fusivities using equations ^ and (O, and the results 
ar e illustrated in Figure [2] Similar to the figures shown 
in IWardld (j2007f ) , we mark different magnetic diffusion 
regimes with different filling color. Six magnetic diffusion 
regimes are considered depending on the relative order 
among 770, rjH and r]A- It is clear that Hall and AD 
becomes more and more important at smaller density 
(surface layer and large disk radii) and larger magnetic 
field strength, as expected {rjH scales as B/ue while rjA 
scales as B^/rie). The addition of grains dramatically 
changes the pattern in the figure at z < 3H, this is be- 
cause grains carry most of the negative charge instead of 
electrons (due to the extremely low ionization rate) and 
the magnetic diffusivity is largely determined by the less 
mobile ions and grains. In this region, different choices 
of chemical reaction networks can make a big difference 
in the resulting magnetic diffusivity pattern. At disk up- 
per layers (z > 3iJ), grains play a less important role on 
the pattern of magnetic diffusivities painted in the Fig- 
ure since free electrons overwhelms the grains, although 
the ionization fraction is still affected by the grains. 

White vertical lines in Figure [5] show the contours of 
constant effective recombination time (equation (|17p ). 
We see that in all of the four cases, the recombination 
time is at least one order of magnitude smaller than 
the dynamical time scale^. This result looks counter- 
intuitive since the recombination time may be expected 
to be smaller in the disk upper layers due to the low 
gas density. However, this is compensated by the en- 
hanced electron abundance near the disk surface (as 
trc ''^ l/nxe)- Together with the extremely low ionization 

^ Note that our definition of ijr^j^ in equation l |17|l captures the 
most rapid recombination process. It is typically shorter than 
the che mical equilibrium time estimated in fPerez-Becker fc Chian3 
ll2011al ). which is sensitive to the slowest chemical processes 
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Fig. 2. — Regimes of non- ideal MHD effects in fiducial PPD models plotted as contours in the plane of disk height and magnetic field 
strength. The left (right) two panels correspond to 1 (10) AU of a MMSN model, and the upper (lower) two panels correspond to the case 
without grains (with 1% of well mixed 0.1/im grain). Different regimes of non-ideal MHD effects arc painted with different background 
colors. Red and magenta: Ohmic resistivity dominated; Dark and light blue: Hall effect dominated; Green and yellow: AD dominated. 
Subdivisions of the color scheme are indicated in the plots. Black contours show constants of the Elsasser number Atot which is increased 
by factors of 10 from bottom left to upper right, with the Atot = 1 contour marked in bold. The bold dashed black line indicates where the 
magnetic pressure equals to the gas pressure (/3 = 1). White vertical lines correspond to contours of constant effective recombination time 
t<=ff . 
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level in PPDs, this result demonstrates that the strong 
coupling limit applies in essentially most regions of typ- 
ical PPDs, and it justifies that single fluid treatment of 
the gas dynamics in PPDs is generally sufficient. In par- 
ticular, the single fluid simulations by BSll on the effect 
of AD on the MRI is d irectly relevant t o PPD s, while 
two-fluid simulations by iHawlev fc Stond ()1998l ) are not 
quite applicable. 

Moreover, we emphasize that our conclusion that 
^f^ch ^ 1 is obtained by using the complex chemical 
network. The usage of a simp le network such as the 
iQppenheimer fc Dalgarnol ()1974|) model can lead to dif- 
ferent conclusions: At 1 AU without grains, we find that 
ircb is about one order of magnitude longer when cal- 
culated with the simple network, and becomes longer 
than the dynamical time at \z\ < H. This is be- 
cause of the lack of recombination channels in the sim- 



ple network, and is relevant to the "revival" of the 
dead zone by turbulent mixing of free electrons from 
the active layer to the midlane, as seen in multi-fluid 
simulations w ith a co-evolying s imple chemical reac- 
tion network (iTurner et al.l 120071: iTurner fc SmmI 120081: 
lllgner &: NelsonI |2008[ ) . However, with a more realis- 
tic chemical network, the reactivation of the dead zone 
by turbulent mixing would appear less likely to oc- 
cur, because the turbulent eddy t ime is comparable 
to the dynamica l time (jPYo mang fc Paoaloizoul 120061 : 
iTurneret al']l2006t ICarbalhdo et al.l 1201 It ) and most free 
electrons would be swallowed by the more rapid recom- 
bination process before being mixed down to the mid- 
plane. Therefore, the density profile of all charged species 
in PPDs should be close to local ionization equilibrium, 
which justifies our adopted criteria in Section [3.71 
In Figure [21 we also plot contours of constant Elsasser 
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?7tot as 



Here we define the Elsasser number based on 



Afot = 



^^^710 



(27) 



where rjtot = x/Vo + ^h + Va^ ^^'^ ^^ measures the im- 
portance of all non-ideal MHD effects as a whole. The 
Elsasser number is in general larger in upper right and 
smaller in lower left, with the Atot = 1 contours plotted 
in bold black lines. Non-ideal MHD terms dominate the 
induction term when Atot < Ij with the properties of 
the MRI significantly modified. In addition, we plot the 
contour of /3 = 1 in bold dashed line. Before applying 
our criteria ([20]) . which will be the subject of the next 
subsection, we note that the bold solid and bold dashed 
lines may serve as rough boundaries enclosing the MRI 
active region. We see that in all cases the dominant non- 
ideal MHD processes in these regions are the Hall effect 
and AD. 

As we discussed in Section 13. 7[ we have ignored the 
Hall effect on the sustainability of the MRI turbulence. 
Here we discuss the potential caveats from this simplifi- 
cation. We are most interested in the boundary between 
MRI active and inactive regions. The lower boundary 
set by A = 1 is relatively well constrained based on the 
study of ISano fc Stone (20023), although more numeri- 
cal study in the Hall dominated regime is still necessary 
since it is generally the case that tjh > Vo near the lower 
boundary. The main uncertainty comes from the upper 
boundary (/? > j3^in{Am)) , which is based on simulations 
of BSll that include only AD. From Figure[2l we see that 
AD is indeed the dominant non-ideal MHD effect close 
to the line of /3 = 1 at 10 AU. Therefore, our criterion 
/3 > /3min(^TO) should provide a relatively reliable upper 
boundary on the strength of the magnetic field for MRI 
to operate. At 1 AU with grains, AD is also likely to be 
the main limiting factor on the magnetic field strength, 
while in the grain-free case, the location of the upper 
boundary is likely to be in the Hall dominated regime, 
and our criterion may require modification. 

4.3. Active Layer and Accretion Rate 

Using the profile of magnetic diffusivities in the pre- 
vious subsection, we estimate the location of the active 
layer by applying our criteria (|20p . In Figure [31 we show 
the A = 1 contour (black bold solid) as well as contours 
of constant Am (black thin solid) in the z-B plane similar 
to Figured! The former (Ohmic resistivity) controls the 
lower boundary of the active layer, while its exact loca- 
tion is determined by the strength of the magnetic field: 
lower for strong field and higher for weak field. This is 
because the Ohmic Elsasser number v\/r]o^ is propor- 
tional to B^, and the disk is better ionized (smaller r/o) 
in the upper layer. The upper limit on the magnetic field 
strength characterized by (3 (blue solid) is controlled by 
AD, where smaller Am value requires weaker magnetic 
field. The value of Am above the lower boundary is gen- 
erally greater than 0.1, and reaches up to 10 in the disk 
upper lay ers, which agrees with PAH-fre e model calcu- 
lations in iPerez-Becker fc Chiangl (|2011a[ ) , and the min- 
imum value of (3 ranges from about 10 to 800. The two 
constraints combined together give the MRI permitted 



region in the z-B plane at a given radius the PPD, as 
painted in gray. 

We see that in all of the four cases, maintaining an 
MRI active layer is possible if the magnetic field is not 
too strong. That is, the required field strength to over- 
come the Ohmic resistivity in the disk upper layer is 
generally much weaker than the equipartition field. This 
is partly due to the fact that rjo oc 1/xe decreases to- 
ward the disk surface more rapidly than density (at least 
for z < AH), making the A = 1 contour steeper than 
the contours of /3 = constant. Also, Am in the disk up- 
per layer is generally greater than 1, where the required 
field reduction from equipartition strength is only mod- 
erate. The MRI permitted region in the z-B plane shifts 
to disk upper layers when sub- micron grains are included 
because of increased r]o ■ It extends toward the disk mid- 
plane as one moves to disk outer regions because the 
reduction of disk surface density slows the recombina- 
tion process and allows ionization photons/particles to 
penetrate deeper. At 10 AU without grains, even the 
disk midplane becomes active for field strength of a few 
times O.OIG. These results are all consistent with previ- 
ous chemistry calculations. 

For MRI-driven accretion, the required magnetic field 
strength at various accretion rates from equation (|25p is 
also shown in Figure [3] (red dashed). We see that in the 
grain-free calculations, the maximum field strength in 
the MRI permitted region corresponds to M « 10~^Mq 
yr~^, while in the presence of sub- micron grains, the 
maximum accretion rate is dramatically reduced to well 
below 1O^^M0 yr^^. The reduction is mostly due to the 
retreat of the lower boundary: the A = 1 contour shifts 
upward, reducing the maximum allowed field strength 
because of reduced gas pressure; in addition. Am be- 
comes smaller, reducing the allowed field strength fur- 
ther. 

We note that most previous calculations either adopt 
the magnetic Reynolds number Rcm = cl/rjo^ = 100 
as the boundary between the active layer and the dead 
zone, which implicitly as s umes /? = 100 in the entire 
disk (iFromang et al.ll2002l: lllgner fc NelsonI 120061: BG09; 
iPerez-Becker fc Chiang! 1201 lal lbl). or adopt the Elsasser 
number criteri on, but assume some constant magnetic 
field strength (!Turner fc Drake! 12009! ). From Figure |3 
the ReAf = 100 criterion corresponds to using the in- 
tersection between the bold solid line (A = 1) and thin 
blue line (/3 = 100) as the boundary. We see that this 
criterion roughly applies in the grain-free case, but may 
substantially overestimate the active column when grains 
are present (i.e., the critical ReM should be larger than 
100). In other words, the reduction of the active layer 
column density by grains is more serious than previously 
estimated, complementing our discussion in the previous 
paragraph. 

Our graphical illustration of the maximum accretion 
rate based on the approximate formula ()25p can be im- 
proved with a more quantitative calculation to give an 
absolute upper limit of the accretion rate. This upper 
limit can be obtained from equation (|23p . where the in- 
tegral is performed across the entire disk height that the 
MRI permitted region extends to, with B chosen to be 
the maximum value determined by /3niin(^"^)- For the 
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Fig. 3. — Similar to Figure[2l but for constraints on the MRI permitted region in PPDs. The bold solid contours correspond to the Ohmic 
Elsasser number A = 1, while the thin solid curves show the contours of constant Am labeled by log jg (Am), terminated at Am = 0.1 . 
Blue bold and thin lines mark the plasma /3 = 1 and /3 = 100 respectively. Permitted regions for the MRI in PPDs based on criteria II20I I 
are painted in gray. The red dashed lines indicate the required field strength in the MRI permitted region corresponding to accretion rate 



of 10 ^,10 ^ and 10 ^ Mq yr ^, respectively, based on equation l|25| l. 

MMSN disk model, wc have 



Mr. 



Cs 



Bt 



dz 



0.98 X lO^^Au'' 



z>0 

active 



B„ 



G 



' dz 1 

^ Mq yr-i . 

(28) 



For the four calculations in Figure [31 wc find the maxi- 
mum accretion rate (in unit of Mq yr~^) to be 4.5 x 10~^, 
1.5 X 10~* for the grain free model at 1 AU and 10 AU 
respectively, 3.8 x 10~^^ and 2.5 x 10^^" for well-mixed 
sub-micron grain model at 1 AU and 10 AU. As expected, 
they agree with the estimate from our graphical illustra- 
tion. 

Our equation (|23p provides an convenient way for 
estimating the accretion rate from the magnetic field 
strength, but it does not predict as a priori the field 
strength in the disk. Here we consider the following 
thought experiment and hypothesis from which we pro- 



pose a way to predict the strength of the magnetic field 
in the active layer hence the accretion rate. 

Suppose the initial magnetic field strength in the disk 
is sufficiently weak (say, /3 = 10^ at the disk midplane). 
According to Figure O the MRI first develops in the 
very upper layer of the disk. The MRI amplifies the 
magnetic field, and consequently, both the upper and 
lower boundaries of the active layer moves toward the 
midplane because of increased field strength^. The MRI 
ultimately saturates into turbulence, and turbulent stir- 
ring maintains ro ughly constant field stren gth across the 
active layer (e.g., lOishi fc Mac Lowll2009l ). We further 
hypothesize that the magnetic field in global disks is 
able to self-organize into a configuration to maximize the 

^ This is somewhat related to the recurrent growth and decay of 
the MRI turbulence as a result of field a mplification and resistive 
damping observed bv lSimon et al.l l)2011h . who adopted a constant 
resistivity profile. In real situations with a rapidly increasing re- 
sistivity toward the disk midplane, the transition may be a more 
smoothed process. 
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rate of angular momentum transport by the MRI turbu- 
lence. While unjustified due to the lack of systematic 
gl obal studies (most g l obal s imulations of the MRI such 
as iFromang fc Nelsoa (J2006{ ) and local stratified simula- 
tions have zero net vertical fiux, which leads to relatively 
weak turbulence), this hypothesis predicts the magnetic 
field strength and accretion rate that lie on the larger 
side than in reality (i.e., an optimistic estimate). 

Based on the discussion above, we predict the accre- 
tion rate as follows. We scan over the strength of the 
magnetic field B in the active layer {B is assumed to be 
constant across the active layer), where for each B, we 
calculate M using equation ([23]). The optimistically pre- 
dicted accretion rate is the maximum M from the scan. 
For the four calculations in Figure |3l the predicted accre- 
tion rates (in unit of Mg yr'^) are 2.7 x 10"^, 7.7 x lO'^ 
for grain free models at 1 AU and 10 AU respectively, and 
2.3 X 10~^^ and 1.3 x 10~^° for well-mixed sub-micron 
grain models at 1 AU and 10 AU. They are about half 
the value of the upper limits listed before. 

While the framework we just described is fully gen- 
eral, the results reported in this paper are subject to a 
few caveats due to simplifications in the adopted disk and 
ionization models that may result in small uncertainties 
up to a factor of a few. These uncertainties are discussed 
below and can be overcome by adopting more realistic 
disk and ionization models that takes into account ra- 
diative transfer, thermodynamics and magnetic pressure 
support to yield more reliable estimate of the accretion 
rate. 

We have assumed constant vertical temperature pro- 
file in our estimate of the accretion rate, while in real- 
ity, the upper layer of the disk is heated by both the 
X-ray photons which are responsible for the io nization 
(|Glassgold et al.l[20qllGorti fc Holleiibach|2004D and the 
MRI itself (BG09. iHirose fc Turneij 120111 ). raising the 
temperature above the midplane temperature by a fac- 
tor of a few. This effect increases the gas pressure, hence 
plasma /3, and allows stronger magnetic field in the active 
layer, which may lead to higher M than our predicted 
value by a factor of a few. 

We have assumed that the gas density profile is Gaus- 
sian in our calculations. In reality, magnetic pressure 
plays more important role in the hydrostatic equilib- 
rium at the disk surface, a nd local isothermal str atified 
shearing-box simulations by iMiller fc Stond ()2000f ) show 
deviation from Gaussian density profile with gas den- 
sity at the disk surface a factor of up to 10 higher than 
our adopted Gaussian model. This may promote accre- 
tion by allowing stronger magnetic field near the surface. 
However, since the ionization rate depends on the col- 
umn gas density from the disk surface, this effect is com- 
pensated by the reduced ionization rate and enhanced 
recombination rate (at given disk height), giving only or- 
der unity corrections to the predicted accretion rate. In 
fact, the predicted accretion rate does not sensitively de- 
pend on the density distribution across disk height, and 
the app roximate calculations by iPerez-Becker fc Ghiand 
(|2011b[ ) that are free from the vertical density distribu- 
tion are generally consistent with ours wit hin order unity. 

Simulations by iMiller fc Stond (|2000( ) also indicates 



non-zero stress in the magnetic corona, but the a value 
is generally one order of magnitude smaller than that in 
the disk. The coronal contribution to the accretion rate 
is about aPchc, with the coronal gas pressure Pc orders 
of magnitude smaller than that in the disk midplane / 
active layer, and the coronal thickness he on the order of 
a few disk scale height (not much thicker than that of the 
active layer). Therefore, for our purpose, accretion in the 
magnetic d ominated disk corona can be saf ely neglected. 

Recently, iPerez-Becker fc Chiang! ()2011bD called for at- 
tention to the role of the FUV ionization, which pro- 
duces much higher ionization fraction than other ion- 
ization sources (making Am > 1000) with penetration 
depth of about 0.01 g cm^^ < J^p^y < 0.1 g cm'^. We 
note that in the MMSN model, 0.1 g cm~^ corresponds 
to z « 3.8H at 1 AU and z « 2.9H at 10 AU, which 
are very upper layers in the disks. Comparing with Fig- 
ure [3] we see that at such heights, the cquipartition field 
strengths in the active layer correspond to M of a few 
times 1Q~^'^Mq yr"-'^ and a few times lO^^Af© yr"-*^ re- 
spectively, which set the upper limit for the MRI accre- 
tion rate driven by FUV ionization. The upper limit is a 
factor of a few below our predicted M in the grain-free 
case, but may well exceed the predicted rate when small 
grains are included. 

In sum, we have provided a general framework for es- 
timating the location of the active layer as well as the 
accretion rate at a given location of PPDs. The most 
important feature is that it explicitly incorporates the 
dependence on the magnetic field strength, with the field 
strength estimated by our physically motivated hypothe- 
sis. Furthermore, the simulation based relation a = 1/2/3 
allows us to provide an optimistic estimate as well as a 
robust upper limit of the accretion rate, without involv- 
ing additional assumptions about the value of a. This 
framework will be extensively used for parameter study 
of the accretion rate in the next section, and can be gen- 
eralized with more realistic criteria and disk models in 
the future. 



5. PARAMETER STUDY 

In this section, we perform a series of chemistry cal- 
culations at different disk radii from lAU to lOOAU and 
explore a number of different model parameters including 
grain size, ionization rate and disk mass. In each series 
of the parameter study, we derive the MRI-driven accre- 
tion rate M as a function of disk radius using the method 
illustrated in the previous section, in parallel with the 
predicted magnetic field strength. The results are shown 
in Figures 2] and [Sj with detailed description and dis- 
cussion given in the subsections that follows. Obviously, 
using a fixed disk model, the predicted M is not neces- 
sarily constant as a function of disk radius, which would 
lead to non-steady accretion if MRI is the only driving 
mechanism. On the other hand, the predicted spatially 
varying M implies that modifications to the adopted disk 
model is necessary for steady state accretion. Detailed 
modeling of the disk structure to match the steady state 
accretion is beyond the scope of this paper, but the trend 
can be readily obtained from our studies. 
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Fig. 4. — The optimistically predicted accretion rate (in unit of Mq 
MMSN disk with cosmic ray and X-ray ionization, model with X-ray 
ionization (lower left) and model with the Desch's disk (lower right). 
(black circles), with well-mixed 1/xm grains (blue triangles) and with 

5.1. The Effect of Grain Size 

We begin by considering our fiducial model: MMSN 
disk with both cosmic ray ionization and X-ray ionization 
[Lx = lO'^^erg s~^). In the upper left panel of Figure U 
we show the predicted M as a function of disk radius. We 
have considered the grain-free case as well as models with 
1/xm and 0.1 /xm well-mixed grains with mass fraction of 
1%. Again, we remind the reader that our predicted 
accretion rate is optimistic. 

In general, the predicted M increases with disk radius 
r in the inner disk, and decreases with r in the outer disk. 
In the former case, the disk midplane is "dead" and ac- 
cretion is layered. Since M ex aSactiveT/fi (where Sactivo 
is the column density of the active layer) , assuming both 
^active and a to be constant, one obtains M (x r for the 
MMSN temperature profile. In reality, as we discussed 
in Sections 13.71 and 14.31 a is controlled by Am (which 
is about 1-10 in the active layer), while the thickness of 
the active layer depends on the mutual effects of Ohmic 



yr~^) as a function of disk radius for our fiducial model (upper left): 
luminosity 100 times higher (upper right), model without cosmic-ray 

In each panel, we show results for calculations in the grain-free case 
O.lfim grain (red squares). 

resistivity and AD, neither assumptions strictly hold. In 
Figure HI we see that the increase of M with disk radius 
in the inner disk is somewhat slower than r. 

When the entire disk becomes active, Ohmic resistivity 
becomes essentially irrelevant in constraining the accre- 
tion rate (see the upper right panel of Figures [2] and [3] at 
r = 10 AU for reference). The accretion rate is mainly 
controlled by the disk surface density, temperature, and 
the magnetic field strength (M oc aST/fi) through AD. 
For MMSN, ET/i7 ex r~^/^, which gives the main source 
of accretion rate reduction. The predicted profile of M 
is also affected by the radial profile of the a parameter 
determined by Am. 

In the grain-free calculation, the entire disk becomes 
active at r > 4 AU in our fiducial model, and the pre- 
dicted M is between 10"^ and 10"*^ M© yr'^ at aU disk 
radii considered. The reduction of the accretion rate 
caused by dust grains is substantial. At 1 AU, M is 
reduced by about two orders of magnitude in the pres- 
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ence of O.l^m grains, while for 1/i grains the reduction is 
more than one order of magnitude. The entire disk be- 
comes active for r > 70 AU and r > 10 AU respectively 
in the above two cases. By and large, the predicted M in- 
creases with disk radius and flattens out at large radii as 
discussed before. The rise in predicted M in the 0.1 ^m 
case toward 100 AU is mainly because the cosmic-ray 
ionization takes over the X-ray ionization (see next sub- 
section). 

We find an interesting fact that at very large disk radii 
(> 40 AU), calculations with l^m grains lead to com- 
parable or slightly higher M than those without grains. 
By carefully checking the equilibrium abundance of all 
chemical species, wc find that in the grain-free calcula- 
tions, ionization level near the disk midplane is deter- 
mined by inetal abundance, while the number density of 
other ions is much smaller (sec Figure [1] for example). In 
the presence of grains, metals are depleted at low tem- 
perature (see equation (26) of BG09) due to adsorption 
onto grain surfaces, and the dominant ion component 
becomes species such as H;^ and HCO"'" (and others) de- 
pending on the gas density. The suppression of metals is 
one of the reasons for the reduction of ionization level, 
but it appears that at very low density and tempera- 
ture, slightly higher ionization level can be achieved due 
to non-metal ions. Nevertheless, we also note that due 
to uncertainties in the reaction rate coefficients in the 
UMIST database, the equilibrium abundance of various 
chemical specie s such as HCO^ can be uncertain up to 
a factor of -- 4 (jVasvunin et al.ll2008f ). 

We see that even in the absence of dust grains, the 
optimistically predicted M is only comparable or below 
10~^Mq yr~^, which is about the media n value of the ob- 
served accretion rate in T-Taur i stars (|Hartmann et al] 
119981 ISiciha-Aguilar et al.l[2005[ ). In the presence of dust 
grains, the accretion rate that can be driven from the 
MRI is reduced dramatically. The situation is particu- 
larly serious for T-Tauri disks, since the gas has to enter 
through the inner disk to accrete. Given the fact that 
our choice of parameters are typical, these results pose 
strong challenge on the effectiveness of the MRI in driv- 
ing rapid accretion in PPDs. Below we further explore 
other effects to see whether higher MRI-driven accretion 
rate can be achieved. 

5.2. The Effect of Ionization Rate 

Because of the uncertainty in the amount of cosmic- 
rays, which may be shielded by the stellar winds, we 
also perform a scries of chemistry calculations without 
cosmic-ray ionization, and the results are shown in the 
lower left panel of Figure 21 We note that besides deeper 
penetration depth, the cosmic-ray ionization rate gener- 
ally exceeds the scattering component of the X-ray ion- 
ization rate beyond about 25 AU, and one might ex- 
pect that cosmic-ray ionization makes significant contri- 
butions to our predicted M in the outer disk. However, 
this does not seem to be the case. 

In the grain-free case, we find that the the transition 
radius larger than which the entire disk becomes active 
shifts from about 4 AU in the fiducial model to about 
10 AU in the absence of cosmic rays. Nevertheless, the 



predicted M does not change by much throughout the 
disk (only slightly reduced). In the presence of grains, 
the transition radius shifts outward to about 25 AU for 
the 1/im case, and is above 100 AU for the O.l^m case. 
Again, our predicted M is also slightly reduced (gener- 
ally within a factor of 2) but still similar to the fiducial 
results. To explain, we note that M largely depends on 
the magnetic field strength in the active layer. Although 
X-rays do not penetrate as deep as cosmic rays to acti- 
vate the disk midplane in the outer disk, the permitted 
magnetic field strength in the active layer turns out to 
be similar to the case with cosmic rays. According to 
equation (|23)) , given similar field strength, A'l is then de- 
termined by the geometric thickness of the active layer, 
which is on the order of H and is not sensitive to whether 
the active layer extends to the disk midplane or not. 

To examine the role of X-ray ionization, we further 
consider a model with Lx = 10"^^ erg s~^, 100 times 
larger than our fiducial X-ray luminosity. Although 
such X-ray luminosity is unusually high for T-Tauri 
stars (|Giidel et al.ll2007[ ). it provides an upper limit on 
the accretion rate that can be driven by X-ray ioniza- 
tion. Alternatively, one might also view this unreal- 
istic choice of high X-ray luminosity as an incorpora- 
tion of other possible but uncertain strong ionization 
sources such as energetic prot ons from protostellar ac- 
tivities (jTurner fc Drake! |2009[ ). We see from the upper 
right panel of Figure S] that the predicted M is much 
higher than our fiducial model. The rise in M is mainly 
due to deeper penetration in the inner disk (< 5 AU) 
and due to reduced AD in the outer disk. Accretion rate 
in grain-free calculations reaches a few times 1O~^M0 
yr~^, which is close to the upper limit of the observed 
accretion rate. Also, comparing the three curves with 
their counterparts in the fiducial plot (upper left) indi- 
cates that increasing X-ray luminosity is more efficient in 
raising accretion rate when small grains are present. In 
sum, it appears that for the MRI to explain the observed 
rate of accretion in PPDs. much stronger ionization rate 
than fiducial is needed. 

We can also compare our results with the predicted 
M by F UV ionization shown bv iPerez-Becker fc Chiang 
(|2011bf) ^. With the FUV penetration depth of < 0.1 
g cm~^, they predict M of a few times 10~^^Mq yr"-'^ 
at 1 AU and it increases roughly linearly with radius 
since accretion only proceeds in the upper surface layer 
with roughly constant surface column density. Compar- 
ing with their Figure 5 we see that FUV ionization is 
likely to be a comparably important ionization source 
as X-rays (while they drive accretion in different lay- 
ers), and in the presence of small grains, FUV ionization 
makes significant contributions to the total MRI-driven 
accretion rate in the inner disk, and could dominate other 
ionization sources in the outer disk. 



^ Tw o order-unity effects may lead IPerez-Becker fc Chiand 
ll2011bl') to overestimating the accretion rate by a factor of a few 
compared with ours: (a) Their equation (6) underestimates the 
density in the disk upper layer by a factor of > 3 in a Gaussian 
density profile, thus promoting higher ionization; (b) Their calcu- 
lations correspond to the maximum possible accretion rate I I28I1 . 
while our predicted rate is about 2 times smaller (see Section 14.311 . 
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5.3. The Effect of Disk Mass 

We explore the role of disk mass by repeating our 
calculations using the Desch's disk model, while the 
adopted ionization rates remain fiducial, and the results 
are shown in the lower right panel of Figure U) We 
see that at 1 AU, where the Desch's disk has about 30 
times more surface density than the MMSN disk, the pre- 
dicted M for all three cases (with and without grains) 
are smaller than the fiducial results. In the grain- free 
case, the predicted M rapidly increases with radius in 
the inner disk, and peaks at about r = 10 AU where 
the entire disk becomes active. The peak rate exceeds 
the MMSN model at the same location by a factor of 
several. The predicted rate falls off towards the outer 
disk more rapidly than that in the MMSN disk, which 
is essentially because the surface density of Desch's disk 
decreases with radius more rapidly. With the grains, the 
predicted M generally increases with disk radius until 
saturation (the entire disk becomes active), and the cor- 
responding accretion rate in the outer disk is comparable 
or slightly exceeds that in the MMSN model. Finally, the 
rates approach the fiducial results at around 100 AU sim- 
ply because the surface densities in the two disk models 
are comparable. 

Our results suggest that given the same ionization 
sources, a more massive disk leads to slower accretion if 
the accretion is layered, while it is able to sustain faster 
accretion if the entire disk is active. The reason is that at 
the same column density from the disk surface (hence the 
same ionization rate) , the gas density is higher in a more 
massive disk (as one can check with the error function) . 
Therefore, the recombination rate is enhanced, reducing 
the ionization fraction hence M roughly as p~^l'^ . On the 
other hand, higher disk mass in the disk permits stronger 
magnetic pressure in the active layer, and dominates the 
former effect when the entire disk column is active. 

We have seen that the predicted M by the MRI tur- 
bulence generally increases with radius in the inner disk, 
and decreases with radius in the outer disk, with the 
transition radius depending on the dust content. This 
result is obviously inconsistent with steady state accre- 
tion. Further evolution would lead to pile-up of mass to- 
wards the inner disk, with the density profile in the outer 
disk becoming more flattened. According to the results 
in this subsection, the pile-up of mass in the inner disk 
would lead to slower rather than faster accretion, which 
further enhances the mass pile-up. Therefore, the fact 
that M decreases with disk surface density in the inner 
disk leads to a runaway pile-up of mass, which may be 
unsustainable. This result further support the conclu- 
sion that mechanisms other than the MRI is needed to 
drive rapid accretion through the inner disk. 

5.4. Predicted Magnetic Field Strength 

In Figure [SJ we show the (optimistically) predicted 
magnetic field strength in terms of the ratio of midplane 
gas pressure to the magnetic pressure in the active layer 
Ps = Pmid/ Pmag- This is a physically more convenient 
quantity than the absolute magnetic field strength and 
provides a guide to numerical modelers. The magnetic 



(29) 



field strength can be obtained by 

B = I3l3;^/\llj^^^ G (MMSN) 
B = 72l3;^/\^lj^^ G (Desch) 

Note the steep dependence of the magnetic field strength 
on disk radius. 

Because the active layer resides in the upper disk, with 
lower gas pressure hence smaller magnetic field, /3s is gen- 
erally much larger than 1. We see that /3s spans a large 
range between about 100 and 10^, depending on the sur- 
face density, ionization rate and disk radius. The pre- 
dicted field strength is generally stronger (smaller /3s) in 
the grain-free case since the active layer resides higher in 
the disk surface in the presence of grains, with lower gas 
pressure hence smaller magnetic field. The field strength 
can be raised significantly by strong ionization (upper 
right panel) , accompanied by increased accretion rate as 
discussed in Section 15.21 For the massive disk model 
(lower right), /3s in the inner disk is much higher than 
the fiducial case because the active layer resides at a 
higher position than in the fiducial model due to the fi- 
nite ionization penetration depth. 

It has been predicted that the presence of ordered mag- 
netic field in PPDs can lead to grain alignment with 
the magnetic field and p roduce polarized emi ssion in 
the millimeter continuum ()Cho fc Lazarianll2007l ). Grain 
alignment is expected in the outer disk, where radiative 
torque dom inates thermal collisions. Recently, however, 
iHughes eta l. (2009) reported non-detection of polarized 
emission using arcsecond-resolution Submillimeter Array 
(SMA) polarimetric observations. One suspected reason 
has to do with the magnetic field strength, which has to 
be above some critical field streng th for the a lignment 
to occur. For 10 — lOO/im grains, iHughes et a l. (200!|) 
calculated the critical field strength to be 10 — lOOmG 
at the location of 50 — 100 AU (corresponding to the an- 
gular resolution of the SMA for the observed sources). 
However, at such distances, our predicted field strength 
is less than 2mG (taking /3s = 100). Therefore, the non- 
detection of polarized emission is, in fact, consistent with 
the presence of the MRI, while if polarized emission were 
detected, the implied magnetic field would be too strong 
for the MRI to operate in the outer disk. 

6. SUMMARY AND DISCUSSION 

We have explored the non-ideal MHD effects in proto- 
planetary disks (PPDs) using chemistry calculations with 
a complex reaction network including both gas-phase and 
grain-phase reactions. Cosmic-ray and X-ray ionization 
processes arc included with standard prescriptions. The 
equilibrium abundance of all charged species are used to 
calculate the full conductivity tensor, from which diffu- 
sion coefficients for Ohmic resistivity. Hall effect and am- 
bipolar diffusion (AD) are evaluated. One major finding 
from our chemistry calculations is that the recombination 
time is much shorter than the orbital time in essentially 
all regions in PPDs, no matter grains are present or not. 
Together with the extremely low level of ionization in 
PPDs, this verifies the applicability of the "strong cou- 
pling" limit, and the gas dynamics in PPDs can be well 
described in a single-fluid framework with magnetic dif- 
fusion coefficients in non-ideal MHD terms given from 
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Fig. 5. — The ratio of midplanc gas pressure to the optimal magnetic pressure at the active layer as a function of disk radius for our 
fiducial model (upper left): MMSN disk with cosmic ray and X-ray ionization, model with X-ray luminosity 100 times higher (upper 
right), model without cosmic-ray ionization (lower left) and model with the Desch's disk (lower right). In each panel, we show results for 
calculations in the grain- free case (black solid), with well-mixed 1/jm grains (blue dashed) and with 0.1/jm grain (red dash-dotted). 

chemical equilibrium. In particular, turbulent mixing of 
chemical species (especially electrons) can be compen- 
sated by the rapid recombination process. Additionally, 
we have updated the calculation of the sticking coeffi- 
cient in electron-grain collisions where grain charge is 
taken into account. The electron sticking coefficient can 
deviate substantially from 1 which has a strong effect on 
the electron recombination rate and grain charge distri- 
bution, and should be treated with care. 

Using the magnetic diffusion coefficients from the 
chemistry calculation, we estimate the location and ex- 
tent of the regions in PPD models where the MRI can op- 
erate to drive disk accretion (i.e., the active layer). Our 
adopted criteria are based on the Ohmic Elsasser number 
A, where A > 1 is required for the active layer as been 
shown by previous sim ulations (e.g.. iTurner et all 120071 : 
lllgner fc Nelsonll2008f ) . and the AD Elsasser number Am, 
where /3 > f3^in{Am) is required for sustained turbu- 
lence from the most recent study by BSl l. We have ig- 
nored the Hall effect based on the study bv lSano fc Stonel 



()2002bD . although the Hall regime is still yet to be more 
carefully explored with numerical simulations. Unlike 
most previous studies, we have considered the depen- 
dence of the diffusion coefficients (translated to the El- 
sasser number) on the magnetic field strength and show 
that the magnetic field strength can be a main limit- 
ing factor on the extent of the active layer because of 
AD. Our study shows that the conventional magnetic 
Reynolds number criterion with Rcm > 100 for the MRI 
to operate may substantially overestimates the column 
density of the MRI active layer in the presence of small 
grains. 

We provide a general framework for estimating the up- 
per limit of the MRI-driven accretion rate. Furthermore, 
we hypothesize that the MRI amplifies the magnetic field 
to maximize the accretion rate, from which we are able 
to make optimistic predictions on the accretion rate as 
well as the magnetic field strength in PPDs. Using this 
framework, we run a series of chemistry calculations with 
different model parameters to study the location and ex- 
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tent of the MRI active layers in PPDs, and to study the 
dependence of MRI-driven accretion rate M on various 
parameters. The main resuhs are summarized as fohows. 

1. Active layer always exists. Its upper boundary is 
set by AD, requiring the field strength not to be 
too strong. The lower boundary is generally deter- 
mined by the Ohmic resistivity and the magnetic 
field strength, with the latter connected to the field 
strength in the upper boundary. 

2. The predicted M increases with disk radius in the 
inner disk where accretion is layered, and flattens 
or decreases with radius in the outer disk when 
the disk midplane becomes active. The transition 
radius depends on grain abundance, and is 4 — 10 
AU in the grain-free case and above 50 AU in the 
presence of well-mixed sub-micron grains. 

3. For our standard model, the predicted M is a few 
times lO~^Af0 yr ~^ in the grain-free case, and 
is reduced by one to two orders of magnitude in 
the presence of sub-micron grains in the inner disk. 
Reduction of M by grains is less significant in the 
outer disk. 

4. The MRI-driven accretion rate is sensitive to the 
protostellar X-ray luminosity but insensitive to 
the deeper-penetrating cosmic-ray ionization rate. 
Extremely strong X-ray luminosity or additional 
strong ionization source with sufficient grain de- 
pletion is needed to achieve accretion rate of '^ 
10-^Mq yr-i. 

5. In the inner disk where accretion is layered, the pre- 
dicted M increases with disk radius, but decreases 
with increasing disk surface density. This situation 
would lead to runaway mass pile-up in the inner 
disk if the MRI were the only mechanism for driv- 
ing accretion in the inner region of PPDs. 

6. The midplane gas pressure is generally a factor of 
100 to 10^ times higher than the predicted mag- 
netic pressure in the active layer. The ratio is 
smaller for stronger ionization and higher for larger 
disk mass or in the presence of small grains. The 
predicted magnetic field strength in the outer disk 
is consistent with the non-detection of polarized 
emission resulting from grain alignment. 

Our results place strong constraints on the effective- 
ness of the MRI-driven accretion in PPDs, especially for 
the inner disk with radius of about 1 — 10 AU where 
accretion is likely to be layered. There are two main dif- 
ficulties concerning points 3 and 5 above: (a) The opti- 
mistically predicted accretion rate under standard model 
prescriptions is about one order of magnitude or more 
smaller than the typical observed value; (b) MRI-driven 
accretion would lead to the runaway pile-up of mass in 
the inner disk. The first difficulty might be alleviated 
by incorporating additional (but uncertain) ionization 
processes such as energetic protons from the protostel- 
lar activities. The second difficulty docs not appear to be 



easily reconciled, since external ionization sources always 
lead to layered accretion. The accumulation of mass may 
trigger the gravitational instability in the early phase of 
PPDs, which may further lead to outburst and ep isodic 
accretion as FU Ori-like events (jZhu et al.l l2009| ) . In 
later phases, mass accumulation may sir nply leads to ver y 
large surface density in the inner disk (|Zhu et al.ll2010l) . 
Since the inner disk (< 10 AU) is smaller than the reso- 
lution of_ttie_cuTreiitly_available sub-millimeter observa- 
tions (jAndrews et al.ll2Q09t) . whether mass accumulation 
occurs may be distinguishable in future observations such 
as ALMA. 

One possible resolution to the above difficulties may 
be ac hieved by a magnetized w ind from the disk sur- 
face (jBlandford fc Pavnd [l98l . The wind could be 
launched from the disk surface where the magnetic field 
is too strong for the MRI to operate, transporting an- 
gular momentum vertically by the Maxwell stress via 
ordered magnetic fields. Non-ideal MHD effects has 
been i ncorporated i nto wind models (jWardle fc Koenigl 
Il993t iTeitleiJ 1201 1[ ). and under favorable conditions, 
wind-driven accretion rate can becom e substantial 
(jKoniglet al.l 120101: ISalmeron et~alll2011t ). It has been 
Toposed that disk w ind may co-exist with the MRI 
Salmeron et al.l [20071 ) . while they operate at different 
vertical locations. Current wind models for PPDs gen- 
erally have a number of unconstrained free parameters, 
and further investigation (especially by numerical simu- 
lations) is needed to assess the effectiveness of disk wind 
on PPD accretions. 

Our results are also applicable in the gas dynam- 
ics of transitional disks, where the MRI is likely to 
be r esponsible for driving accret ion from the outer 
disk (jChiang fc Murrav-Clavl[2007h . with the inner gap 
/ hole opened by m ultiple planets that gu i de the 
gas streams th rough (|Perez-Becker fc Chiand l2011al : 
IZhu et all 1201 It ) . Our predicted accretion rate in the 
outer disk with r > 10 AU is on the order of 1O~^M0 
yr~^, with relatively weak dependence on the dust con- 
tent compared with the inner disk. The result is 
consistent with th e range of observed accretion rate 
(jNaiita et al.ll2007l ) , thus MRI alone appears sufficient to 
account for the accretion rate in transitional disks. This 
conclusion can als o be achieved with far UV ion ization 
instead of X-rays (jPerez-Becker fc Chiangir2011b[ ) . 

Our study of the MRI active layer in PPDs represents 
a first attempt to estimate the MRI-driven accretion rate 
that is based on the most up-to-date non-ideal MHD sim- 
ulations. In the mean time, it is limited by the predictive 
power of the simulations themselves. In particular, the 
non-linear properties of the MRI in the Hall dominated 
regime is still unexplored. Moreover, the AD simulations 
by BSll are unstratified, and whether their obtained cri- 
terion for sustained MRI turbulence holds remains to be 
verified in stratified simulations. In conclusion, our study 
calls for further progress in non-ideal MHD simulations, 
and in turn, more realistic criteria for sustained MRI 
turbulence can be easily incorporated into our general 
framework after future simulations. 
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APPENDIX 
STICKING COEFFICIENT IN ELECTRON-GRAIN COLLISIONS AND GRAIN CHARGING 

The collision of an electron with a grain particle does not necessarily lead to absorption because the electron has 
excess energy as it reaches the grain surface. The absorption probability P^ depends on the size of the grain a, grain 
charge Z, and the electron energy E. The sticking coefficient is defined as the absorption probability averaged over a 
thermal distribution of electron energies, which is a function of temperature T 

where <t{E) = a{a,Z,E) is the electron-grain cross section is given bv iDraine &: Sutml (|1987f ) (see their Section II). 
The task here is to calculate Pe{E) = Pe{a, Z,E). 

The stick i ng pr obability in the case of electron collision with neutral grains was derived in the Appendix B of 
iNishi et al.l (|1991[ ). Below we generalize their derivation to include grain charge. 

Let Eq be the initial energy of the electron (at infinity), D be the depth of the potential well between the electron 
and the grain due to the electric polarization interaction (on the order of a few eV) . The kinetic energy of the electron 
as it reaches the grain surface is therefore Eq + D + Ze'^/a. In an inelastic collision with a surface atom/molecule 
(whose mass is Ms) of the grain, the upper limit of energy transfer is given by 

. „ 4?7ip , „ Ze^ „, , . , 

AE.u^—^{Eo + +D) , A2 

Ms a 

where rrie is the electron mass. From the phonon theory, the probability of inelastic collision electron collision a, and 
the averaged energy transfer in each inelastic collision SE, are given by ([Umebavashi &: Nakanol[l980t ) 

where To is the Debye temperature of the lattice and is about 420K for graphite. The electron has to experience 
about I inelastic collisions on average before being truly absorbed, where / is given by 

ISE > Eo>il- 1)SE . (A4) 

Between (say, the ith and the (i -I- l)th) inelastic collisions, the electrons undergo elastic collisions with the gra in and 
has a probability /3i to escape. Therefore, the absorption probability is given by (|Umebayashi fc Nakanol 119801 ) 



P.= 



'' 1-A 



Swrrfe^- ^^^) 



The remaining task is to calculate Pi, corresponding to the electron escape probability when its total energy is 
Ei ~ Eq — iSE. This is made convenient by considering energy conservation to write the radial energy of the electron 
as 

E, = E, + + -^^^ ^ - - {E, + + D) sin^ , (A6) 

where the 2nd and 3rd terms represent the potential energy due to net charge and electric polarization, and the last 
term represents the transverse energy obtained from angular momentum conservation L^/2mr^, with 9 being the angle 
between the velocity vector immediately after the elastic collision (scattering) and the radial direction. We see that 
Er approaches infinity close to the grain surface due to electric polarization, and approaches Ei at infinity. If Er 
becomes negative at some intermediate radius, it means that the electron can not travel beyond this radius and has 
to return to the grain surface for another collision. If E^ is positive at any r, the electron trajectory is open and will 
escape. Therefore, the purpose is to find the critical value of 0cr, at which the function Er{r) touches the horizontal 
axis tangentially: Er = and dE^/dr ~ 0. The latter condition can be more conveniently replaced by d{r'^Er)/dr = 0. 
To proceed, we introduce dimensionless variables x = r/a, A = Da/e^ and e^ = Eia/e^, and obtain 

'^ + 2.2 fal _ n + ! = ^(^'^ + ^ + ^) ^-' '- ' (A^) 
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Fig. 6. — Electron sticking coefficient for O.lfim (left) and 1/im (right) sized grains as a function of gas temperature. Two values of the 
binding energy D=leV (solid lines) and D=3eV (dashed lines) are chosen. In each case, we show curves for different grain charges with 
Z = —4 (blue), Z = (black), Z = 2 (green), and Z = A (red). Their sticking coefficients monotonically decrease as Z becomes larger and 
more positive. We adopt D=leV for all our chemistry calculations. 

2xe, + Z^ J" . (A8) 

One can solve for x i n equatfon (|A8p. which always has one solution for a; > 1. We note that equation (|A8p is equivalent 
to equation (2.10) of iDraine fc SutinI ()1987() in the calculation of cross section, which is natural since electron escape 
is the inverse process of electron capture. Back substituting the solution of x to equation (jA7[) one obtains the critical 
angle 6'cr- 

One caveat in the calculation is that the left hand side of equation (jATj can become negative when Z is negative and 
ej is sn iall^. In fact, this expression is proportional to the electron-grain cross section (equation (2.8) of lDraine fc SutinI 
(|l987t )). This situation correspond to that a potential well with Er < always exist which prevents the electron from 
approaching the grain (in the case of electron capture) , or escaping from the grain (in the case of we are interested 
here). If this occurs, we simply set the escaping probability /3i = (or 6'cr = 0). In fact, this means that after the ith 
collision, the electron is destined to be adsorbed even if its total energy Ei > 0. 

Finally, by considering the process of elastic collisions as isotropic scattering, one obtains the absorption coefficient 
to be 



^^=2^ 



27r sin 



= 1 — cos ( 



(A9) 



As an example, we show the electron sticking coefficient for collisions with O.l^m and 1/im grains in Figure |6l We 
see that s^ monotonically decrease with disk temperature and drops by about 2 orders of magnitude from lOK to 
lOOOK. This is because the electron has more excess energy at higher temperature and has to undergo more inelastic 
and elastic collisions. Higher binding energy D raises the sticking coefficient, by reducing the critical angle 9ci at 
each elastic collision. Our new result is that the sticking coefRcient for negatively charged grains is in fact larger than 
that for positively charge grains^. This is reasonable because for negatively charged grains, the electron has been 
decelerated as it approaches the grain, thus bounces fewer times before being absorbed. Finally, the sticking coefficient 
is larger for larger grains. Meanwhile, the difference in the sticking probability for grains with different charged states 
is less for larger grains. Both are understandable since the electric potential due to grain charge at the grain surface 
reduced for large grains. 
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